# conda install -c conda-forge osmnx

# https://osmnx.readthedocs.io/en/stable/
# https://qiita.com/moshi/items/e383491664d028cd2166
# 
# https://zenn.dev/homing/articles/f9a314841c737d
# https://qiita.com/naozo-se/items/bc168a2282907a6b29bf


# https://stackoverflow.com/questions/62285134/how-to-fill-water-bodies-with-osmnx-in-python

import osmnx as ox
import networkx as nx
import folium
from geopy.distance import distance

# 与えられた緯度経度

latitude = 34.69155489576708  # 例としてケーニッヒスベルクの緯度
longitude = 135.48959938567188  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
G = ox.graph_from_point((latitude, longitude),
                        dist=distance_km*1000)
   #network_type='drive')

# グラフを取得
nodes, edges = ox.graph_to_gdfs(G)

# 橋の情報を抽出（'bridge'タグを持つエッジのみを抽出）
#bridge_edges = edges[edges['bridge'] == 'yes']
bridge_edges = edges

# グラフのノードとエッジを抽出
bridge_graph = nx.Graph()

# エッジ（橋）をグラフに追加
for idx, data in bridge_edges.iterrows():
    u, v, key = idx
    bridge_graph.add_edge(u, v)

# オイラー閉路の判定
eulerian = nx.is_eulerian(bridge_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# グラフのエッジ（橋）を地図に追加
for u, v in bridge_graph.edges():
    point_u = (G.nodes[u]['y'], G.nodes[u]['x'])
    point_v = (G.nodes[v]['y'], G.nodes[v]['x'])
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook

import osmnx as ox
import networkx as nx
import folium

# 与えられた緯度経度
latitude = 54.7088  # 例としてケーニッヒスベルクの緯度
longitude = 20.5108  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
G = ox.graph_from_point((latitude, longitude), dist=distance_km*1000, network_type='all')

# 無向グラフに変換
G = G.to_undirected()

# グラフの連結成分を取得
connected_components = list(nx.connected_components(G))

# 橋の情報を抽出（'bridge'タグを持つエッジのみを抽出）
bridge_edges = [edge for edge in G.edges(data=True) if 'bridge' in edge[2] and edge[2]['bridge'] == 'yes']

# 各連結成分を一つのノード（陸地）と見なす
land_nodes = {i: component for i, component in enumerate(connected_components)}

# 橋のノードを含む陸地ノードのIDを見つける関数
def find_land_node(node):
    for land_id, nodes in land_nodes.items():
        if node in nodes:
            return land_id
    return None

# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()

# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)

# 橋エッジを追加（橋が異なる陸地を接続している場合）
for u, v, data in bridge_edges:
    land_u = find_land_node(u)
    land_v = find_land_node(v)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v, **data)

# オイラー閉路の判定
eulerian = nx.is_eulerian(final_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# 橋のエッジを地図に追加
for u, v, data in final_graph.edges(data=True):
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook

print(land_graph)
print(final_graph)

MultiGraph with 5819 nodes and 8690 edges
Graph with 1 nodes and 0 edges

import folium
import osmnx as ox

import geopandas as gpd

# 対象地域の道路情報取得 (愛知県名古屋市中村区)
#query = "Nakamuraku,Nagoya,Aichi,Japan"
#G = ox.graph_from_place(query, network_type="drive")

(lat,lon)=(34.69199639822432, 135.49084584150586)

# Bounding box
bN,bS,bE,bW=lon+1,lon-1,lon-1,lon+1

#graph=ox.graph.graph_from_point((lat, lon),
#                                dist=1000,
#                                network_type='drive')
graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
                                tags={"natural": "water"})

# 道路グラフネットワーク可視化
#fmap=ox.plot_graph_folium(graph)
fmap=ox.plot_graph_folium(graph)
folium.LayerControl().add_to(fmap)
fmap


#fmap.save(outfile="road_network.html")
#opts = {"node_size": 5, "bgcolor": "white", "node_color": "blue", "edge_color": "blue"}
#ox.plot_graph(G, show=False, save=True, filepath="road_network.png", **opts)

/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/1645265506.py:18: FutureWarning: The `geometries` module and `geometries_from_X` functions have been renamed the `features` module and `features_from_X` functions. Use these instead. The `geometries` module and function names are deprecated and will be removed in the v2.0.0 release. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
/opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:48: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  return features.features_from_bbox(north, south, east, west, tags=tags)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/array.py:950, in GeometryArray.estimate_utm_crs(self, datum_name)
    949 try:
--> 950     return CRS.from_epsg(utm_crs_list[0].code)
    951 except IndexError:

IndexError: list index out of range

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Cell In[17], line 18
     13 bN,bS,bE,bW=lon+1,lon-1,lon-1,lon+1
     15 #graph=ox.graph.graph_from_point((lat, lon),
     16 #                                dist=1000,
     17 #                                network_type='drive')
---> 18 graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
     19                                 tags={"natural": "water"})
     21 # 道路グラフネットワーク可視化
     22 #fmap=ox.plot_graph_folium(graph)
     23 fmap=ox.plot_graph_folium(graph)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:48, in geometries_from_bbox(north, south, east, west, tags)
     22 """
     23 Do not use: deprecated.
     24 
   (...)
     45 gdf : geopandas.GeoDataFrame
     46 """
     47 warn(DEP_MSG, FutureWarning, stacklevel=2)
---> 48 return features.features_from_bbox(north, south, east, west, tags=tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:134, in features_from_bbox(north, south, east, west, bbox, tags)
    131 polygon = utils_geo.bbox_to_poly(bbox=bbox)
    133 # create GeoDataFrame of features within this polygon
--> 134 return features_from_polygon(polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:349, in features_from_polygon(polygon, tags)
    346 response_jsons = _overpass._download_overpass_features(polygon, tags)
    348 # create GeoDataFrame from the downloaded data
--> 349 return _create_gdf(response_jsons, polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:440, in _create_gdf(response_jsons, polygon, tags)
    437 relation_types = {"boundary", "multipolygon"}
    439 # extract geometries from the downloaded osm data
--> 440 for response_json in response_jsons:
    441     response_count += 1
    443     # Parses the JSON of OSM nodes, ways and (multipolygon) relations
    444     # to dictionaries of coordinates, Shapely Points, LineStrings,
    445     # Polygons and MultiPolygons

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/_overpass.py:389, in _download_overpass_features(polygon, tags)
    373 """
    374 Retrieve OSM features within boundary from the Overpass API.
    375 
   (...)
    386     a generator of JSON responses from the Overpass server
    387 """
    388 # subdivide query polygon to get list of sub-divided polygon coord strings
--> 389 polygon_coord_strs = _make_overpass_polygon_coord_strs(polygon)
    390 utils.log(f"Requesting data from API in {len(polygon_coord_strs)} request(s)")
    392 # pass exterior coordinates of each polygon in list to API, one at a time

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/_overpass.py:253, in _make_overpass_polygon_coord_strs(polygon)
    234 """
    235 Subdivide query polygon and return list of coordinate strings.
    236 
   (...)
    249     list of strings of exterior coordinates of polygon(s)
    250 """
    251 # first subdivide the polygon if its area exceeds max size
    252 # this results in a multipolygon of 1+ constituent polygons
--> 253 poly_proj, crs_proj = projection.project_geometry(polygon)
    254 multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
    255 multi_poly, _ = projection.project_geometry(multi_poly_proj, crs=crs_proj, to_latlong=True)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/projection.py:60, in project_geometry(geometry, crs, to_crs, to_latlong)
     57     crs = settings.default_crs
     59 gdf = gpd.GeoDataFrame(geometry=[geometry], crs=crs)
---> 60 gdf_proj = project_gdf(gdf, to_crs=to_crs, to_latlong=to_latlong)
     61 geometry_proj = gdf_proj["geometry"].iloc[0]
     62 return geometry_proj, gdf_proj.crs

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/projection.py:99, in project_gdf(gdf, to_crs, to_latlong)
     97 # else if to_crs is None, project gdf to an appropriate UTM zone
     98 elif to_crs is None:
---> 99     to_crs = gdf.estimate_utm_crs()
    101 # project the gdf
    102 gdf_proj = gdf.to_crs(to_crs)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/geodataframe.py:1451, in GeoDataFrame.estimate_utm_crs(self, datum_name)
   1415 def estimate_utm_crs(self, datum_name="WGS 84"):
   1416     """Returns the estimated UTM CRS based on the bounds of the dataset.
   1417 
   1418     .. versionadded:: 0.9
   (...)
   1449     - Prime Meridian: Greenwich
   1450     """
-> 1451     return self.geometry.estimate_utm_crs(datum_name=datum_name)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/geoseries.py:1214, in GeoSeries.estimate_utm_crs(self, datum_name)
   1178 def estimate_utm_crs(self, datum_name: str = "WGS 84") -> CRS:
   1179     """Returns the estimated UTM CRS based on the bounds of the dataset.
   1180 
   1181     .. versionadded:: 0.9
   (...)
   1212     - Prime Meridian: Greenwich
   1213     """
-> 1214     return self.values.estimate_utm_crs(datum_name)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/array.py:952, in GeometryArray.estimate_utm_crs(self, datum_name)
    950     return CRS.from_epsg(utm_crs_list[0].code)
    951 except IndexError:
--> 952     raise RuntimeError("Unable to determine UTM CRS")

RuntimeError: Unable to determine UTM CRS

import osmnx as ox
ox.settings.log_console = True

# add more items here to get all the landforms you want
places = ['Manhattan, NY, USA', 'Brooklyn, NY, USA', 'Queens, NY, USA', 'Bronx, NY, USA']
land = ox.geocode_to_gdf(places)

# get the water bodies
left, bottom, right, top = land.total_bounds
bbox = top, bottom, right, left
poly = ox.utils_geo.bbox_to_poly(*bbox)
water = ox.features_from_polygon(poly, tags={'natural': 'water'})

# constrain the plotting window as desired
c = land.unary_union.centroid
bbox = ox.utils_geo.bbox_from_point((c.y, c.x), dist=12000)

water_color = 'blue'
land_color = '#aaaaaa'
fig, ax = ox.plot_footprints(water, bbox=bbox,
                             color=water_color, bgcolor=water_color,
                             show=False, close=False)
ax = land.plot(ax=ax, zorder=0, fc=land_color)

/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/273932149.py:11: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  poly = ox.utils_geo.bbox_to_poly(*bbox)

import osmnx as ox
import folium

# 場所を設定（例：東京都渋谷区）
place_name = "Shibuya, Tokyo, Japan"

(lat,lon)=(34.69199639822432, 135.49084584150586)

# 陸地と道路のグラフを取得
#G = ox.graph_from_place(place_name, network_type='drive')
G = ox.graph_from_point((lat,lon), network_type='drive')

nodes, edges = ox.graph_to_gdfs(G)

# 中心点を取得
lat = nodes['y'].mean()
lon = nodes['x'].mean()

# foliumマップを作成
m = folium.Map(location=[lat, lon], zoom_start=14)

# 陸地を描画
folium.GeoJson(data=edges).add_to(m)

# マップを保存または表示
m.save('map.html')

# マップを表示（Jupyter Notebookの場合）
m

Make this Notebook Trusted to load map: File -> Trust Notebook

import osmnx as ox
import folium
import geopandas as gpd

# 場所を設定（例：東京都渋谷区）
place_name = "Shibuya, Tokyo, Japan"

# 陸地のポリゴンを取得
land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})

# 水域のポリゴンを取得
water_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': ['water', 'bay', 'strait']})

# 道路のグラフを取得
G = ox.graph_from_place(place_name, network_type='drive')
nodes, edges = ox.graph_to_gdfs(G)

# 中心点を取得
lat = nodes['y'].mean()
lon = nodes['x'].mean()

# foliumマップを作成
m = folium.Map(location=[lat, lon], zoom_start=14)

# 陸地を描画
for _, polygon in land_polygons.iterrows():
    folium.GeoJson(polygon['geometry'], style_function=lambda x: {'fillColor': 'green', 'color': 'green'}).add_to(m)

# 水域を描画
for _, polygon in water_polygons.iterrows():
    folium.GeoJson(polygon['geometry'], style_function=lambda x: {'fillColor': 'blue', 'color': 'blue'}).add_to(m)

# 道路を描画
folium.GeoJson(data=edges).add_to(m)

# マップを保存または表示
m.save('map.html')

# マップを表示（Jupyter Notebookの場合）
m

/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/388463500.py:9: FutureWarning: The `geometries` module and `geometries_from_X` functions have been renamed the `features` module and `features_from_X` functions. Use these instead. The `geometries` module and function names are deprecated and will be removed in the v2.0.0 release. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})

---------------------------------------------------------------------------
InsufficientResponseError                 Traceback (most recent call last)
Cell In[22], line 9
      6 place_name = "Shibuya, Tokyo, Japan"
      8 # 陸地のポリゴンを取得
----> 9 land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})
     11 # 水域のポリゴンを取得
     12 water_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': ['water', 'bay', 'strait']})

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:128, in geometries_from_place(query, tags, which_result, buffer_dist)
    104 """
    105 Do not use: deprecated.
    106 
   (...)
    125 gdf : geopandas.GeoDataFrame
    126 """
    127 warn(DEP_MSG, FutureWarning, stacklevel=2)
--> 128 return features.features_from_place(query, tags, which_result, buffer_dist)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:298, in features_from_place(query, tags, which_result, buffer_dist)
    295 utils.log("Constructed place geometry polygon(s) to query API")
    297 # create GeoDataFrame using this polygon(s) geometry
--> 298 return features_from_polygon(polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:349, in features_from_polygon(polygon, tags)
    346 response_jsons = _overpass._download_overpass_features(polygon, tags)
    348 # create GeoDataFrame from the downloaded data
--> 349 return _create_gdf(response_jsons, polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:489, in _create_gdf(response_jsons, polygon, tags)
    487 if len(geometries) == 0:  # pragma: no cover
    488     msg = "No data elements in server response. Check log and query location/tags."
--> 489     raise InsufficientResponseError(msg)
    491 # remove untagged elements from the final dict of geometries
    492 utils.log(f"{len(geometries)} geometries created in the dict")

InsufficientResponseError: No data elements in server response. Check log and query location/tags.

import osmnx as ox
import networkx as nx
import folium

# 与えられた緯度経度
latitude = 54.7088  # 例としてケーニッヒスベルクの緯度
longitude = 20.5108  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
G = ox.graph_from_point((latitude, longitude), dist=distance_km*1000, network_type='all')

# 無向グラフに変換
G = G.to_undirected()

# 橋の情報を抽出（'bridge'タグを持つエッジのみを抽出）
bridge_edges = [(u, v) for u, v, data in G.edges(data=True) if 'bridge' in data and data['bridge'] == 'yes']

# 橋のノードを収集
bridge_nodes = set()
for u, v in bridge_edges:
    bridge_nodes.add(u)
    bridge_nodes.add(v)

# 橋を除いたグラフを作成
land_graph = G.copy()
land_graph.remove_edges_from(bridge_edges)

# 陸地の連結成分を取得
land_components = list(nx.connected_components(land_graph))

# 各連結成分を一つのノード（陸地）と見なす
land_nodes = {i: component for i, component in enumerate(land_components)}

# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()

# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)

# 橋エッジを追加（橋が異なる陸地を接続している場合）
for u, v in bridge_edges:
    land_u = next((land_id for land_id, nodes in land_nodes.items() if u in nodes), None)
    land_v = next((land_id for land_id, nodes in land_nodes.items() if v in nodes), None)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v)

# オイラー閉路の判定
eulerian = nx.is_eulerian(final_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# 橋のエッジを地図に追加
for u, v in final_graph.edges():
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook

import networkx as nx
import matplotlib.pyplot as plt

# final_graphの可視化
plt.figure(figsize=(10, 8))
pos = nx.spring_layout(final_graph)  # ノードの位置を計算

# ノードを描画
nx.draw_networkx_nodes(final_graph, pos, node_size=300, node_color='lightblue')

# エッジを描画
nx.draw_networkx_edges(final_graph, pos, width=2, edge_color='blue')

# ノードのラベルを描画
nx.draw_networkx_labels(final_graph, pos, font_size=12, font_family='sans-serif')

# グラフを表示
plt.title("Final Graph Visualization")
plt.axis('off')
plt.show()